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Abstract. We consider a chain composed of N coupled harmonic os- 
cillators in contact with heat baths at temperature Tg and T r at sites 1 
and N respectively. The oscillators are also subjected to non-momentum 
conserving bulk stochastic noises. These make the heat conductivity sat- 
isfy Fourier's law. Here we describe some new results about the hydro- 
dynamical equations for typical macroscopic energy and displacement 
profiles, as well as their fluctuations and large deviations, in two simple 
models of this type. 



The time evolution and the nature of the resulting non-equilibrium sta- 
tionary states (NESS) of systems in contact, at their boundaries, with ther- 
mal reservoirs at different temperatures remains a challenging problem. The 
only such systems with bulk Hamiltonian dynamics for which NESS are 
known are harmonic crystals [TJI2]- These are described by specifying the 
positions and momenta of the oscillators, qj and pj, j = 1, 2, . . . , N. The 
Hamiltonian for a D = 1 chain (with obvious extension to D > 1) is given 



where we have set the mass of each particle and the nearest neighbor coupling 
equal to 1. We impose boundary conditions qo = qN+i = 0, po = Pn+i = 0. 
The strength of the pinning potential is regulated by the parameter v > 0. 
This system is put in contact with thermal reservoirs at temperatures T# at 
site 1 and T r at site N via the Ornstein-Uhlenbeck (Langevin) process. 

In the absence of any stochastic bulk noise the heat flux in the non- 
equilibrium stationary state is independent of N, corresponding to the heat 
conductivity diverging linearly in N p] . To imitate the effects of interactions 
between phonons and to produce NESS with normal transport one may add 
bulk noise to the system [3]-[8]. The new NESS will, we expect, exhibit 
certain universal features also present in realistic anharmonic systems. It is 
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the purpose of this note to summarize some new results about such systems. 
A more detailed analysis will be given elsewhere [6]. 

We shall now consider two ways of adding noise to the system. The 
first case we consider is the so-called self-consistent model [BJ. Each site 
j G {1, . . . , N} is connected to a Langevin reservoir at temperature Tj > 0. 
The temperatures of the reservoirs at the boundary sites are fixed by the 
conditions T\ = Ti, Tjv = T r , while the temperatures of the interior site 
reservoirs are determined self-consistently by requiring that in the NESS 
there is no net flux of energy between the system and the interior reservoirs. 

The second model we consider is one in which we add in the bulk flip 
dynamics which consists of reversing the velocity of each particle at random 
independent Poissonian times. These flips are energy but not momentum 
conserving and we call the model the velocity flip model @]. 

Let ju(q, p, t) be the probability density of our system at time t. The time 
evolution of fx is described by the Fokker-Planck equation 
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where B\. ^ represents the bulk noise with b 
b = s for the self-consistent case: 



v for the velocity-flip and 
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, k = 2,...,N- 1. 



The NESS of the two models, denoted by u s>Y and u s>s respectively, are 
solutions of (jl.2p with, for the self-consistent model, 7\, k = 2, . . . , N— 1, 
determined by the self-consistency condition /i SjS (p|) = (p|) s ,s = T k . 

If the temperatures Ti and T r are equal to a common value T, the steady 
state of the self-consistent chain and the steady state of the velocity flip 
model are both equal to the Gibbs state with temperature T that we denote 
by (•^i' 71 ). This is a Gaussian measure with covariance C eq (T) = TC eq (l). 
For Ti ^ T r , the NESS of the two models are different. However, it is easy to 
see that both stationary states are centered, and by the results derived in [4], 
we know that the two point correlation functions of both models coincide, 
when the rate of flipping and the coupling to the internal reservoirs are the 
same. Nevertheless, the two measures are very different: /x S)S is Gaussian but 
/x S)V is a non-trivial mixture of Gaussian states. This has a simple non-trivial 
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consequence: If /(q, p) is a linear function of the form 

N 

/(q,p) = ^2(aiqi + biPi) , 

i=l 

where (a±, . . . , ajsr, b\, . . . , 6jv) are arbitrary real numbers, then 

(/ 4 )s,v>3[(/ 2 ) s , v ] 2 = (/ 4 ) s , s . 

This follows from the fact that for a centered Gaussian measure (/ 4 ) = 
3 [(/ 2 )] , and (-) s ,v is a superposition of such measures. 

Both models satisfy Fourier's law with the same value of the conductivity. 
Indeed, since the microscopic energy current across the bond (k — l,k) 

3k = ~2^ Pk +Pk-\){<lk ~ Qk-i), k = 2,. . . ,N - 1 , 

is a quadratic function of the momenta and positions, its averages over (-) S;S 
and (-)s,v are equal. The macroscopic current J = lim7v->oo [kn(T£ — T r )/N], 
with k = lim at_>oo kjv , is given according to [3] by 

1/7 

1.3 



2 + v 2 + 0^2 + 4) 

The outline of the rest of the paper is as follows: In Section [2] we present 
the macroscopic equations for the density profiles of the bulk conserved 
quantities in the velocity-flip model. In Section [3] we show that long-range 
energy correlations are present in the velocity-flip model but not in the self- 
consistent one. Section [4] contains numerical simulations of these long-range 
correlations. 

2. HYDRODYNAMICAL SCALING LIMIT OF THE VELOCITY FLIP MODEL 

We have to distinguish two cases according to whether v = (unpinned) 
or v > (pinned). The unpinned case is similar to that investigated by one 
of the authors in [5] . 

2.1. The unpinned chain. When v = the bulk dynamics conserves two 
quantities. The first one is the energy T~L. The second one is the deformation, 
r x of the lattice, with r x = q x +i — q x , x = 0,1, ... ,N. This second 
conservation law has to be taken into account in the hydrodynamic analysis. 
The energy at a site x E {1, . . . , TV} is now given by 

222 

fy] iy* ! y* 

and £q = \r 2 , £n+i = \ r N- To establish the hydrodynamic limits corre- 
sponding to the two conservation laws, we look at the process with time 
scaled by iV 2 and space scaled by N, i.e., in the diffusive scale [8j. Assume 
that initially the process is started with a Gibbs local equilibrium measure fi 
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associated with a macroscopic deformation profile uo(q) and a macroscopic 
energy profile £o(q): 

. N+l 

A = ^ II exp {-^(x/iV)^ - r (x/iV)r x .)} , (2.2) 

where To = Pq 1 and ro are the temperature and tension profiles correspond- 
ing to the given energy and deformation profiles assumed to be continuous. 
Then we have for any macroscopic point q £ [0, 1] 

Jim (r\ Nq ](0)) = u (q), lim (£\ Nq ](0)) = e Q (q) , (2.3) 

where [y] is the integer part of y and the averages are w.r.t. /t. 

The question then is: what happens at any later (macroscopic) time i? 
It is shown in [6] that 

lim (r [Nq] (N 2 t)) = u(q, t), lim (£ [Nq] (N 2 t)) = e(q, t) , (2.4) 
where u, e are solutions of the following macroscopic diffusion equation 



dtu = 7 1 d 2 u 
3 t e = (27)" 1 d\ (e + u 2 /2) 



-l & ( B _l ..2 ( 2 - 5 ) 



with the initial conditions u(q,0) = uo(q), e(q,0) = £o(q). Eq (|2.5p is to be 
solved subject to the boundary conditions 

d q ll(0,t) = d q V,(l,t) = 0, 

e--Uo,t)=T i , ( E --Ul,t) = T r . 

The boundary condition on d q u comes from the fact that the length of the 
chain remains fixed. 

Taking the limit t — > oo in these equations we obtain the typical macro- 
scopic profiles of the system in the NESS, i.e., a flat deformation profile 
u = and a linear profile T interpolating between Ti and T r , 

e(q)=T(q)=T e + (T r -T e )q, (2.7) 

for the energy profile. 

2.2. The pinned chain. Assume that the system is initially distributed 
according to a Gibbs local equilibrium measure associated to the energy 
profile £o(q), q 6 [0, 1], and define e(q, t) as the evolved profile in the diffusive 
scale, i.e., 

e(q,t) = lim (£ [Nq] (tN 2 )) . 
Then e is the solution of the following heat equation 

d t e = d q (nd q e) , 

e(q,0) = e (q) , (2.8) 
e(0,t)=T e , e(l,t)=T r . 
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The conductivity k, which is independent of the temperature, is given by 

When t goes to infinity e(q, t) converges to the linear profile T(q) (given in 
(|2.7p ) both for the velocity- flip and the self-consistent model. We note finally 
that, since the self-consistent model does not conserve energy in the bulk, 
we do not expect any autonomous macroscopic equations in that model. 

3. Energy Fluctuations 

3.1. The pinned chain. Our goal is to estimate the probability that in 
the stationary state the empirical energy profile, N (q), defined by looking 
at the microscopic energy S x , for x equal to the integer part of Nq, is close 
to a prescribed macroscopic energy profile e(q) different from T(q), i.e., we 
want to find the large deviation function (LDF) for the NESS. 

At equilibrium = T r = T = the stationary state (-) s . s coincides 
with the usual Gibbs equilibrium measure /i^' eq and by the usual large 
deviations theory (see e.g. [TT]) we have that for any given macroscopic 
energy profile e(-) 

4 m (0 N (<1) ~ e(q)) ~ e-^M (3.1) 
where the large deviation function (LDF) 

V cq (e) = [ 
Jo 

coincides with the difference between the free energy of the system in local 
thermal equilibrium (LTE) and the true equilibrium free energy with e(q) = 
T. 

Out of equilibrium (Tg ^ T r ) there is also a large deviation principle [9] 

<v (*"<•) ~ ~ e-^W 

but the explicit form of V is in general unknown. What is true however, is 
that V depends only on two macroscopic quantities: the heat conductivity 
and the mobility [9]. 

For the system we are interested in, the conductivity k(T) is given by 
(jl.3p and is independent of T. By the Einstein relation the mobility x(T) is 
equal to x(T) = k(T)u(T) where cr(T) is the static compressibility defined 
by the equilibrium correlation 

a(T) = ]r<(£ - T)(S X - T))^ T \ (3.2) 

A simple computation shows that for our system cr(T) = T 2 and Theorem 
6.5 of [9] applies. It follows that V(-) is given by 
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where F is the unique increasing solution of 

( d 2 F F-e 

\jd^f~^^' (3.4) 
(F(0)=T e , F(l) = T r . 

Surprisingly, the function V is independent of the pinning value v 2 and 
of the intensity of the noise 7. In fact, it coincides with the LDF of the 
Kipnis-Marchioro-Presutti (KMP) model considered in [10]. In that model 
the dynamics are entirely stochastic. 

It is now easy to derive the Gaussian fluctuations of the empirical energy. 
We consider a small perturbation, e = T + Sh, of the stationary profile e. 
The functional V has a minimum at T so that 

V{e) = V{f) + l -b 2 (h, C~ l h) + o{5 2 ) 

The operator C is the covariance for the Gaussian fluctuations of the em- 
pirical energy under the invariant measure fi^ v . The computations are the 
same as in [10] and we get 

C = T 2 l + (T r -T e ) 2 (-A )- 1 (3.5) 

where Ao denotes the Laplacian with Dirichlet boundary conditions on [0, 1]. 

3.2. Unpinned Chain. For the unpinned case, we are not able to obtain 
the expression of the LDF of the two conserved quantities nor for the energy 
alone. We conjecture that the LDF for the energy is the same for the pinned 
and unpinned case. Moreover, we are able (under suitable assumptions) to 
show that, in the unpinned case, (|3.5p is still valid [6]. 

3.3. Self-Consistent Model. We show in [6] that the energy correlations 
in the NESS of a self-consistent chain are not long range in the sense that the 
variance of the total energy divided by N is given by the LTE measure. The 
remainder is (at least) of order N^ 1 / 2 in the pinned case and iV- 1 / 4 In 2 N 
in the unpinned case. Energy fluctuations are thus dominated by the local 
equilibrium term. Moreover, in the pinned chain the NESS is reached expo- 
nentially fast on a microscopic time-scale. Therefore, the energy profile is 
then given by the linear equilibrium profile for all macroscopic times. 

4. Numerical simulations 

We have performed numerical simulations in the velocity-flip model, both 
with and without bulk pinning. In these simulations, the initial position and 
momentum of each particle was chosen randomly, from a centered uniform 
distribution in the intervals]— y/2T/k, y / 2T/k] and [— v2Tm, V2Tm] where 
T = r '~^ Tr and k is the interparticle harmonic potential, and we have set 
m = 1 and ks = 1 (the mass and Boltzmann's constant respectively). The 
numerical solution to the dynamics was carried out by choosing a small step 
St, and at each step updating the position and momentum of the particles 
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100 
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0.1 


1.40 


0.01 
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8,1 


0.1 


1.39 


0.01 


1.20 


1.40 


400 
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0.1 


1.42 


0.02 







Table 1. Total energy variation in the unpinned model (y = 0) 
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200 


8,1 


0.25 


0.1 
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200 


8,1 


0.25 


1.0 


1.39 
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1.20 


1.40 
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8,1 


0.25 


0.1 


1.39 


0.01 
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8,1 


0.25 


0.1 


1.46 


0.05 






200 


5,1 


0.25 


0.1 


1.30 


0.01 


1.15 


1.30 



Table 2. Total energy variation in the pinned model {v > 0) 



first according to the Hamiltonian evolution, then for the heat baths, and fi- 
nally taking into account the velocity flips. The solution to the Hamiltonian 
part was obtained by performing an Euler integration, implemented through 
the velocity- Verlet algorithm. A numerical approximation was used to de- 
termine the white-noise forces for the Langevin heat baths, and a standard 
scheme was implemented to generate the Poisson processes for the velocity 
flips at every site. The state of the system was judged to be sufficiently 
close to the steady state when the currents were essentially constant 
throughout the chain. 

The total energy T~L is given by (jl.ip . and our goal is to estimate numeri- 
cally its total fluctuations, by measuring the observable 



s N = N 



(H) 



2 

s.v 



Using equation (|3.5p the measured fluctuations in the total energy should 
be given by 
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local equilibrium approximation, and Coo = ^(T r — X^) 2 / T(q)dq 

correction due to the long-range correlations. The prefactor ^ is o 
by integrating over ((— Ao)~ 1 l)(q) - 
explicit integrals yields then 
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To get good statistics, several realizations of the initial data were con- 
sidered. The results from numerical simulations for the variation in total 
energy, denoted by sjy, are collected in Tables CD and [21 These results are 
compared with the theoretical estimate, Sqo, and with s|^ c ' eq '. We see that 
there is a very close match between the predicted Sqo and the measured 
values. 



The NESS for the self-consistent crystal with coordinates ... s i D ,i a = 
1,2, . . . , N , with thermal reservoirs with temperatures Ti (T r ) in contact 
with all oscillators having values i± = 1 (i\ = N) and periodic boundary 
condition in the transverse directions was solved in [3]. As shown in [4] 
this has the same covariances as the corresponding velocity flip model. In 
particular the heat conductivity k is given by 



where P(y) 2 = is 2 + 2 X^i=i 1 ( 1 - cos(2vn/i)). The same argument which 
gives fj, B v as a superposition of Gaussians in ID also holds in D > 1. The 
hydrodynamical equation for the pinned velocity flip case will again be that 
of the form (I2T5D 



where n is given by (|5.1|) . Unfortunately we do not have any information 
about the LDF in the D > 1 case. 
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